function R = compute_residual_image_median(image, k)
if nargin < 2
    k = 1;
end
R = double(image);
for i = 1:k
    R = residual_image(R);
end
end

function R = residual_image(image)
    p_m = medfilt2(image);
    % after padding, edge values are distorted.
    R = image(2:end-1,2:end-1) - p_m(2:end-1,2:end-1);
end